{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example of surface traction integration for FEM (Bathe book, p. 361)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "using SymPy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\left[ \\begin{array}{r}0\\\\0.333333333333333\\\\0\\\\0.333333333333333\\\\0\\\\1.33333333333333\\end{array}\\right]$\n"
      ],
      "text/plain": [
       "6×1 Matrix{Sym}:\n",
       "                 0\n",
       " 0.333333333333333\n",
       "                 0\n",
       " 0.333333333333333\n",
       "                 0\n",
       "  1.33333333333333"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Analytical integration using SymPy\n",
    "r, p1, p2 = symbols(\"r, p1, p2\")\n",
    "N =  [r*(1+r) 0; 0 r*(1+r); -r*(1-r) 0; 0 -r*(1-r); 2*(1-r^2) 0; 0 2*(1-r^2)]\n",
    "t =  [0; (1+ r)*p1+(1-r)*p2] \n",
    "int_Nt = integrate.(1/4*N*t, r)\n",
    "int_Nt.subs(r,1).subs((p1,p2).=> (1,1)) .- int_Nt.subs(r,-1).subs((p1,p2).=> (1,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5555555555555556"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Define integration points coordinates and weights\n",
    "nip =  3\n",
    "ipx = zeros(nip,1)\n",
    "ipw = zeros(nip,1)\n",
    "# Location\n",
    "ipx[1] = -.775 #+.5\n",
    "ipx[2] = 0.0# +.5\n",
    "ipx[3] = .775 #+.5\n",
    "# Weight    \n",
    "ipw[1] = 5.0/9.0\n",
    "ipw[2] = 8.0/9.0\n",
    "ipw[3] = 5.0/9.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Surface nodes shape functions and derivatives\n",
    "nnfa = 3 # assume quadratic element\n",
    "N    = zeros(nip,nnfa,1)\n",
    "dNdx = zeros(nip,nnfa,1)\n",
    "for i=1:nip\n",
    "    ξ = ipx[i]\n",
    "    if nipfa==3\n",
    "        N[i,:,:]    .= [ξ/2*(ξ-1)  1.0-ξ^2 ξ/2*(ξ+1)]'\n",
    "        dNdx[i,:,:] .= [ξ-1.0/2.0 -2.0*ξ   ξ+1.0/2.0]'  #w.r.t η2\n",
    "    end\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6-element Vector{Float64}:\n",
       " 0.0\n",
       " 0.0\n",
       " 0.0\n",
       " 0.3336805555555555\n",
       " 1.3326388888888885\n",
       " 0.33368055555555554"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Traction vector\n",
    "nipfa = nnfa\n",
    "T = zeros(3,2)\n",
    "T[:,2] .= 1.0 # only in y\n",
    "Fex     = zeros(8)\n",
    "Fey     = zeros(8)\n",
    "Fs      = zeros(2)\n",
    "Ff      = zeros(2*nnfa)\n",
    "f2n     = [2 5 1]\n",
    "xn      = [-1.0 0.0 1.0]\n",
    "yn      = [-1.0 1.0 1.0]\n",
    "x       = [xn; yn]'\n",
    "for ip=1:nipfa \n",
    "    Ni    = N[ip,:]\n",
    "    dNdXi = dNdx[ip,:,1] \n",
    "    J     = x' * dNdXi\n",
    "    detJ  = 1.0#sqrt( J[1]^2 + J[2]^2 )\n",
    "    \n",
    "    # Surface tractions \n",
    "    Fs[1] = Ni'*T[:,1]  # use shape functions to evaluate traction at integration point\n",
    "    Fs[2] = Ni'*T[:,2]\n",
    "    # Integrate surface traction contributions\n",
    "    Ff[1:nnfa]       .+= ipw[ip] .* detJ .* Ni .* Fs[1]\n",
    "    Ff[(nnfa+1):end] .+= ipw[ip] .* detJ .* Ni .* Fs[2]  \n",
    "end\n",
    "Fex[f2n] .= Ff[1:nnfa]'\n",
    "Fey[f2n] .= Ff[(nnfa+1):end]'\n",
    "# display([Fex; Fey])\n",
    "display(Ff)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.2",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.2"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
